# Full text of "The Hamiltonian structure and Euler-Poincaré formulation of the Vlasov-Maxwell and gyrokinetic systems"

## See other formats

The Hamiltonian structure and Euler-Poincare formulation of the Vlasov-Maxwell and gyrokinetic systems J. Squire, 1 H. Qin, 1 - 2 W. M. Tang, 1 and C. Chandre 3 ^Plasma Physics Laboratory, Princeton University, Princeton, New Jersey 08543, USA 2) Dept. of Modern Physics, University of Science and Technology of China, Hefei, Anhui 230026, China ^Centre de Physique Theorique, CNRS - Aix-Marseille Universite, Campus de Luminy, 13009 Marseille, France We present a new variational principle for the gyrokinetic system, similar to the Maxwell- Vlasov action presented in Ref. Ill The variational principle is in the Eulerian frame and based on constrained variations of the phase space fluid velocity and particle distribution function. Using a Legendre transform, we explicitly derive the field theoretic Hamiltonian structure of the system. This is carried out with a modified Dirac theory of constraints, which is used to construct meaningful brackets from those obtained directly from Euler- Poincare theory. Possible applications of these formulations include continuum geometric integration techniques, large-eddy simulation models and Casimir type stability methods. 1 I. INTRODUCTION An inherent difficulty in studying the dynamics of magnetized plasmas is the enormous separa- tion of important time-scales present in many physical systems of interest. Nonlinear gyrokinetic theory has become an indispensable tool in these inquiries, as it removes the fastest time-scales from the system, while keeping much of important physics relevant to turbulent transport". A particularly nice way to construct a gyrokinetic theory, pioneered in Refs. |5KZL is to use Lie- transforms to asymptotically change into co-ordinates in which gyro-orbit dynamics are decoupled from the rest of the system. A great advantage of this technique, aside from the entirely systematic and formal procedure, is that the single particle equations are guaranteed to be Hamiltonian, with associated conservation properties. Going further, it is advantageous from both a philosophical and practical standpoint to derive the entire system, including both electromagnetic fields and par- ticles, from a single field-theoretic variational principle. These ideas were explored by Sugama 8 and Brizard 9 , who derived gyrokinetic action principles starting from Maxwell- Vlasov theories, as well as in previous work in Refs. 10-13 Some advantages of this type of formulation are a much simplified derivation of the gyrokinetic Maxwell's equations and exact energy-momentum conservation laws through Noether's theorem. Field theories often admit many different varia- 12 17|) . each with its own advantages and tional principles (e.g., for Maxwell- Vlasov see Refs. disadvantages. A good example is the difference between Lagrangian and Eulerian actions; the former being constructed in variables that follow particle motion and the latter in variables at fixed points in phase space. It is interesting to explore new types of variational principles, both for the general understanding of the structure of the theory in question, and for practical applications that may require an action of a particular form. In this work, we present a new gyrokinetic action principle in Eulerian co-ordinates, using Euler-Poincare reduction theory— — on the Lagrangian action in Refs. [s| and 21. In addition, using the reduced Legendre transform and a modified version of the Dirac theory of constraints, we de- rive field theoretic Poisson brackets, similar to the Vlasov -Maxwell 22 26 and Vlasov-Poisson 2 ^ 2 ^ brackets. To our knowledge this is the first explicit demonstration of the field theoretic Hamil- tonian structure of the gyrokinetic system (see Refs. 23 and 28j for a different approach that has recently been used to write down a Poisson brackets for simplified gyrokinetic systems). Our derivation proceeds from the action principle in Ref. Is] and its geometric formulation 2 ^-. We do not purport to derive a gyrokinetic co-ordinate system, but rather formulate the theory based on a given single particle Lagrangian. In this way, it is trivial to extend concepts to deal with more complex gyrokinetic theories, for instance theories with self consistent, time-evolving background fields 21 . We then use the ideas in Ref. 1 to reduce the Lagrangian action to one in Eulerian co-ordinates, based on symmetry under the particle-relabeling map from Lagrangian to Eulerian variables. The variations for this new action principle in the Eulerian frame are constrained, and lead to the Euler-Poincare equations, which are shown to give the standard gyrokinetic Vlasov equation. In some ways the action principle is similar to that of Brizard 9 , in that constrained variations must be used, with both theories having a similar form for the variation of the distribution function, F. Nevertheless, there are significant differences, particularly that our principle is formulated in terms of the Eulerian phase space fluid velocity and is in standard 6-D phase space, rather than 8-D extended phase space. The Eulerian gyrokinetic action of Refs. |29| is quite different to that presented here, with unconstrained variations in 12-dimensional extended phase space and the use of Hamilton-Jacobi functions in the action functional. For more information on this approach to 33 gyrokientic theory, see Refs. 1l?[ [i3, 3oL 3 dre transform is performed 1 , leading straightforwardly to a Poisson bracket. However, this bracket Equipped with the Eulerian action, a reduced Legen- must be reduced to a constraint submanifold before a meaningful form can be obtained, a process that is performed with a modified version of the Dirac theory of constraints^—. Finally, we show how to include the electromagnetic fields in the bracket via second application of Dirac theory. One of our primary motivations in this work is the possibility of utilizing recent ideas from fluid mechanics to develop advanced numerical tools for gyrokinetics. Of particular importance is the idea of geometric integrators, which are designed to numerically conserve various impor- tant geometrical properties of the physical system. For instance, having a numerical algorithm that has Hamiltonian structure can be very important, with profound consequences for the long- time conservation properties 37 . The theory of finite dimensional geometric integrators is relatively well developed 37 , including an application to single particle guiding center dynamics^^. How- ever, many aspects of the construction of field-theoretic geometrical integrators are not as well understood, both for practical implementation and the deeper mathematical theory. One approach, which has yielded fruitful results, is to discretize a variational principle and perform variations on the discrete action to derive an integration scheme. Some examples of field theoretic inte- grators constructed in this way are those for elastomechanics^ 2 -, electromagnetism 43 , fluids and magnetohydrodynamics 4445 , and a particle-in-cell (PIC) scheme for the Vlasov-Maxwell system 4 ^. The results presented in this work would be used to construct a continuum Eulerian gyrokinetic integrator, since our variational principle is in Eulerian form. Analogously, a variational princi- ple in Lagrangian form is used to construct a Lagrangian (particle-in-cell) integrator 46 . We note that in discretizing a variational principle it is obviously not desirable to be in an extended phase space, unless these extra dimensions can somehow be removed after a discretization. As well as integrators, other potential applications of the formulation presented here are the use in stability calculations with Casimir invariants 47 ' 48 and the construction of regularized models for large-eddy simulation 4 ^—. It is significant to note at this point that as the power of modern supercomputing systems contin- ues to advance at a rapid pace toward the exascale (10 18 floating point operations per second) and beyond, it is quite clear that the associated software development challenges are also increasingly formidable. On the emerging architectures memory and data motion will be serious bottlenecks as the required low-power consumption requirements lead to systems with significant restrictions on available memory and communications bandwidth. Consequently, it will be the case in multiple application domains that it will become necessary to re-visit key algorithms and solvers - with the likelihood that new capabilities will be demanded in order to keep up with the dramatic archi- tectural changes that accompany the impressive increases in compute power. The key challenge here is to develop new methods to effectively utilize such dramatically increased parallel comput- ing power. Algorithms designed using geometric ideas could be very important as simulations of more complex systems are extended to longer times and enlarged spatial domains with high physics fidelity. The article is organized as follows. In Section [XT] we clarify the differences between Eulerian and Lagrangian action principles for kinetic theories and explain the Euler-Poincare formulation of the Maxwell- Vlasov system 1 . This is done with as little reference to the formal mathematics as possible, with the hope that readers unfamiliar with the concepts of Lie groups and algebras will understand the general structure of the theory. Section [III] explains the construction of the gyrokinetic variational principle, starting from a given single particle gyrokinetic Lagrangian. We give a brief derivation of the Euler-Poincare equations and show how these lead to a standard form of the gyrokinetic equations. The Hamiltonian structure is dealt with in Section|IV] After formally constructing a Poisson bracket from the Lagrangian, we describe how the modified Dirac theory of constraints is used to reduce the bracket to a meaningful form. Finally, numerical applications are briefly discussed in Section |V] and conclusions given in Section |VT1 Throughout this article we use cgs units. In integrals and derivatives, z denotes all phase space 4 variables, while jc denotes just position space variables. Species labels are left out for clarity and implied on the variables F (or /), m, e, U and M, respectively the distribution function, particle mass, particle charge, Eulerian fluid velocity and momenta conjugate to U. Summation notation is utilized where applicable, with capital indices spanning 1 — > 6 and lower case indices 1 — > 3. II. EULERIAN AND LAGRANGIAN KINETIC VARIATIONAL PRINCIPLES When formulating a variational principle for a continuum fluid-type theory, it is very important to specify whether Lagrangian or Eulerian variables are being used. These notions can be con- fusing in kinetic plasma theories, since one must consider the motion of the phase-space fluid. In addition, unlike the Euler fluid equations, the equations of motion for kinetic plasma theories have the same form in Eulerian and Lagrangian co-ordinates. Considering the Vlasov-Maxwell system for simplicity, a Lagrangian description gives the equation of motion at the position of a particle carried along by the flow (simply a physical particle). One formulates a variational principle in terms of the fields jc (jc , v , t), v (x , v , t), which are the current position and velocity of an ele- ment of phase space that was initially at (xq, Vo). The distribution function is of course just carried along by the Lagrangian co-ordinates, i.e., f(x (x , v , t) , v (jc , v , t)) = f (jc , v ). This type of formulation is the most natural for a kinetic theory, since it is the logical continuum generalization of the action principle for a collection of particles interacting with an electromagnetic field. An Eulerian variational principle is formulated in terms of the velocity of the phase space fluid at a fixed point, U, without the notion of where phase space density has been in the past. Thus, at a point (jc, v), the jc component of the fluid velocity is simply v (the co-ordinate), while the v component is E + v x B/c. The distribution function /, is advected by U, meaning it is the solution to the differential equation d t f = --Cuf = -U ■ V/, where U and V are in six-dimensional phase space. An illustration of these concepts is given in Figure \T\ for the 1-D Poisson-Vlasov system (in 2-D phase space). Finally, we note that in discussing the distinction between Eulerian and Lagrangian actions, we refer only to the plasma component of the variational principle; the electromagnetic fields are always in Eulerian co-ordinates. 5 V n *■ * {x ,v ) // ^ X {x,v) tjj(x ,v ) = U (x,v) FIG. 1. Illustration of the particle relabelling map, i]/(xq, Vq) and its inverse for the one-dimensional Vlasov- Poisson system. A. Euler-Poincare reduction This section gives a very informal introduction to Euler-Poincare theory through a brief review of the Vlasov-Maxwell formulation presented in Ref. 1. We purposefully do not use the group theoretical notation of Refs. fl, [lsi andfl^l (e.g., the ❖ and ad* operations) so as to introduce the general ideas to readers not familiar with Lie groups and algebras. Euler-Poincare type variational principles first appeared in Ref. |20| in the context of magnetohydrodynamics. The purpose of the Euler-Poincare framework is to provide a straightforward method to pass from a Lagrangian to an Eulerian action principle. The important idea is that the La- grangian is invariant under the right action of the particle relabeling transformation, ^(jcq, Vq) = (jc(jco,Vo),v(jc ,Vo)), which maps plasma particles with initial position (jc , v ) to their current position (jc, v). In essence, this invariance means that we can eliminate the extraneous par- ticle labeling information and still have an equivalent system. The map iff acts on / on the right as / = foiff~ l , where f Q is the initial distribution function. This equation is simply / (x (xo,Vo,t) ,v (xo,Vo,t)) = /o(jco, Vo), as discussed above. In the electromagnetic part of the Lagrangian iff does not act on the potentials (p and A, since electromagnetic dynamics must be independent of particle labeling. The phase space fluid velocity in the Lagrangian frame is simply iff(xo, Vq), since this is the rate of change of (jc, v) at the position (jc, v). In contrast, the Eulerian phase space fluid velocity is if/if/~ l , since this operation first takes (jc,v) back to (jc ,v ) with if/' 1 and then gives the velocity at (jc, v) with if/ (jco, Vo), see Figured] The starting point for the reduction is a Lagrangian Lagrangian for the Vlasov-Maxwell system. 6 For instance, L = ^ j dxodvo fo ^- A (jc) + mv j • jc - ^ m ^ 2 - ecf> (jc) + S~n' dX „ 2 -V,-- -|VxA| 2 (1) which is very similar to the action principle of Low 16 . The Vlasov equation follows from the standard Euler-Lagrange equations for \b = (x, v), d 6L 8L _ dt 5\jj Sip (2) along with / (x, v, t) = f (xq, Vq). Maxwell's equations come from the Euler-Lagrange equations for A and (p. This Lagrangian is invariant under the relabeling transformation iff, i.e., L /o ((/r, ip, (p, (p, A, A) = L /or i (lfnfr l ,ifnjr l ,(p, tp,A, A) = l(u,(f>,i>,A,A,F), (3) where U = ij/t/f' 1 is the Eulerian fluid velocity, a vector field. In recognizing that the distribution function is actually a phase space density, we denote F = fdx A dv. Treating F as 6-form rather than a scalar changes the form of certain geometrical operations in the Euler-Poincare equations [Eqs. © and ©] and is very important for the gyrokinetic Euler-Poincare treatment (see Sec- tion Hill). Practically speaking, to construct the reduced Lagrangian, I, one simply replaces (jc, v) with U, and considers jc and v to be co-ordinates rather than fields. Thus, / = ^ J F ^- A + mv j • U x - ^-mv 2 - ecf> (jc) + 8^ ' dx „ dA |Vx A| (4) where U x denotes the jc components of U. The equations of motion are derived from the reduced Lagrangian /, by considering how the unconstrained variations of ib (used to derive the standard Euler-Lagrange equations) translate into constrained variations of U and F. This leads to varia- tions of the form 6U=^-[U,rj], SF = -LqF, (5) dt 7 where rj is in the same space as U and vanishes at the endpoints; and [ , ] is the standard Lie bracket, U.Vj] - tj.VU. Evolution of F is given by the advection equation — +£ V F = 0, (6) ot which arises from the equation / (x, v, t) = / (jco, Vo)- Variation of J dtl with 6U and 6F leads to the Euler-Poincare equations, d 61 „ 61 ^61 where 61/6U is a 1-form density. We give straightforward derivations of Eqs. © and © in Sec- tion UTTA] below. Since F is a 6-form, Xj/i 7 = V • (FU) and Eq. © is the conservative form of the Vlasov equation (see Section ITlI Al for more information). The equations for A and cf> are just the standard Euler-Lagrange equations. Calculation of Eq. © with the Vlasov-Maxwell reduced Lagrangian [Eq. ©] leads to U x = v, U v = E + -v x B, (8) c as expected. The fact that there is no need to solve differential equations for components of U is related to the strong degeneracy in the system (see Sections Hill and HVl below) . To obtain the Hamiltonian or Lie-Poisson form of the equations, one performs a reduced Leg- endre transform as, f 61 h = (M,U)+ dxA- — -l, (9) J 6A where M = 61/6U and the inner product ( , ) is integration over phase space [see Eq. (pTb l. (A thorough treatment of the degeneracies of the system is given in Ref . 1 .) It is then straightforward to show that {r,0} iP = -^ j dzM- 6T_ 6G 6~M , 6~M + f , 3® ^SY 6Y 6S > / dzF\ V V — / \6M 6F 6M 6F . [ j (6Y 6® 6G 6Y\ -4nc / dx\ , (10) J \6A 6E 6A 6E) is an infinite dimensional Poisson bracket for the system; that is, M = {M,h}, F = {F,h}, E = [E,h] and A = {A,h\ are formally the same as the Euler-Poincare equations (using the 8 generalized Legendre transform of Ref. and the Jacobi identity is satisfied. Nevertheless, this manifestation of the bracket has major problems. In particular, the meaning of functional deriva- tives with respect to the M variables can be unclear, since these are constrained due to the linearity of the Lagrangian in U. To overcome these problems and formulate a meaningful bracket on the space of plasma densities and electromagnetic fields, we use a modified version of the Dirac the- ory of constraint s 13134153 . For the convenience of the reader a very brief overview of standard Dirac theory in Appx. |Aj while the modified version used in parts of this work is covered in Appx. |B] The relevant constraints are 0, =M t - F i^-At + mv\ = 0, i = 1 3, 0; = Mj = 0, j = 4 -> 6. (11) We form the constraint matrix Cu (z, z') = {0/ (z) , 0/ (z')} and construct the inverse according to the procedure detailed in Appxs. [A]and|B] This comes out to be, Cjj(z,z') 1 mF (z) 6{z-z')6 ss x -10 1 -1 -^B 7 1 —B z mc * 1 mc > -B x mc A 0-1 ^B Y -^B x mc y mc A (12) which is simply the single particle Poisson matrix (multiplied by 6 (z - z') IF) and B is a function of z. We will see a similar connection to the single particle Poisson bracket in the reduction of the gyrokinetic bracket. We then use Eq. (|A31 ) and restrict the functional Y and to not depend M , i.e., ST/6M = (see Appx.®. The final result, including a change of variables from A to B, is the Poisson bracket for the 9 Mawell-Vlasov system, cm 2 s + 47rV — dT F d@ F d&F BYt. dx dv dx dv dT F dO F dzFB ■ — x — - dv dv + Anc dz dT_ dE K se) 6B d@ ~ dE d r F —F— dv 6B 6® 8E (13) This bracket is identical to that published previously, initially proposed in Ref . 25 with a correction given later in Ref. |22] The derivation above explicitly shows the link between this and the work of Ref. Ill There is a slight taint on the validity of this bracket in that it requires V • B = for the Jacobi identity to be satisfied (as shown in Refs. 23 and 27 through direct calculation). Recently, this obstruction has been partially fixed using the Dirac theory of constraints 24 . Somewhat more detail is given for the derivation of the gyrokinetic bracket (see Section ITVl). which proceeds in a very similar manner. Euler-Poincare reduction is perhaps more natural when applied to fluid systems 19,20 . In this case, there are fewer degeneracies and the Lagrangian/Eulerian distinction is more obviously rel- evant (e.g. one does not measure phase space fluid velocities, in contrast to the Eulerian fluid velocity of a fluid system). Recently, similar ideas have been applied to general reduced fluid and hybrid models in plasma physics^—. III. GYROKINETIC VARIATIONAL PRINCIPLE 21 The Our starting point is the geometric approach to gyrokinetic theory advocated in Ref. general idea is to construct a field theory, including electromagnetic potentials, from the particle Poincare-Cartan 1-form, y. This approach is conceptually very simple; once the interaction of quasi-particles with the electromagnetic field is specified, particle and field equations follow in straightforward and transparent way via the Euler-Lagrange equations. With any desired approxi- mation (e.g., expansion in gyroradius), energetically self-consistent equations are easily obtained without necessitating the use of the pullback operator. The use of these ideas in gyrokinetic simu- lation has been advocated in, for instance Refs. 57 and 58 In this article we consider the particle 10 1-form y as given, its derivation can be found in Refs. |2L|8Und|2H among other works. The Poincare-Cartan form y in 7-D phase space, P, (including time) defines particle motion through Hamilton's equation, i T dy = 0, (14) which is derived from stationarity of the action Jl sp = f y. Here r is a vector field whose integrals define particle trajectories (including the time component) and i denotes the inner product. Note that y is essentially just L sp dt, where L sp is the standard Lagrangian; that is, for y = y a dz a - Hdt, the Lagrangian is simply L sp = y a z a - H. To construct a field theory, y is used to define the Liouville 6-form, Q T = -1-dy Ady Ady. (15) The Liouville theorem of phase space volume conservation is then simply, £, t £It = 0. Introducing the distribution function of particles in phase space /, the field theory action for the interaction of a field of particles with the electromagnetic field is W = 47rJfQ T Ay + J dx£ EM , (16) where -Csm is the electromagnetic Lagrangian density. In this action y is in the Lagrangian frame. For example, in Cartesian position and velocity space, fQ, is in (jco, Vo) co-ordinates, while y is in (jc (jc , v ) , v (jto, v )) co-ordinates. Unless general relativity is important, one can choose r to be of the form r = d/dt + r^, where r z has no time component, and consider 6-D phase space. Defining Q. to be the dX A dP component of Q. T (i.e., no Adt), the dXAdP component of X r ^r is the standard Liouville theorem of phase space volume conservation, -a + £ Tz n = o. (17) We can simplify the variational principle by considering y to be L sp dt and carrying out the wedge product. This type of procedure provides a generalization of the original variational principle of Low^ [e.g., Eq. (OQ)] to arbitrary particle-field interaction. 11 We now specialize to a general gyrokinetic form for the particle Lagrangian, 7 = -A + (X) ■ dX + — udO - Hdt, (18) c e with (19) Here, X is the gyrocenter position, u the gyrocenter parallel velocity co-ordinate, /i the con- served magnetic moment and 6 the gyrophase. The vector field A {X) is the vector potential of the background magnetic field and b {X) is the background magnetic field unit vector. These fields will not be considered variables in the field theory action. The vector R (X) = Vei ■ e%, where e\ (X) J. e 2 {X) J. b (X), is necessary for gyrogauge invariance of the Lagrangian, i.e., invariance with respect to a change in the definition of the 6 co-ordinate. Eq. (TT~8T ) is accurate to first order in e B , the ratio of the gyroradius to the scale length of the magnetic field 2 . The single particle Hamil- tonian, H = jtnu 2 + fiB(X) + H gy , contains both the guiding center contribution, |mw 2 + ^B(X), and the gyrocenter contribution from the fluctuating fields, H gy . For most of this article Hgy will be taken to be a general function of (X, u,fj). Different forms exist in the literature, depending on desired accuracy and fluctuation model used. For instance, in Ref. la H gy is given to second order in eg (the ratio of the magnitudes of the fluctuating fields, 0i and A\, to the background field) as where iff = (p\ — -v ■ A\, p is difference between the particle and gyrocenter positions, and <) denotes an average over 6. The tilde in ^jSi,^}^ denotes the gyrophase dependent part of a function and S i is a gauge function associated with the first order gyrocenter perturbation, (\f/{X + p)>. Eq. (fT8l) is the standard gyrokinetic single particle Lagrangian in Hamiltonian form 2 , meaning all the fluctuating field perturbations are in the Hamiltonian part (dt component) of y. This is the form most suitable for computer simulation 2 ^ and also has the advantage of having the same Poisson structure as the guiding center equations. We now construct a field theory from the single particle Lagrangian Eq. (fT8l . a process often referred to as lifting 23 ^ 60 . We first calculate the phase space component (dX A du A djx A d6 (20) 12 component) of the volume element Q = -\dy A dy A dy. This is simply B^/m = b ■ fm, with ^ = Vx A f , i.e., the standard guiding center Jacobian. In co-ordinates, the variational principle Eq. (PT61) is then simply, ft = / dtU = ^ j dt j dX Q A duo A dfio A dOo—B^fo x|^-A t -X+— + / dtL EM , (21) which is essentially the original gyrokinetic variational principle of Sugama 8 . L EM should be chosen as LeM = ~LJ ^(l V ^l 2 -l Vx ( A + A i)l 2 )' ( 22 ) where x = X + p, so the Ampere-Poisson system is obtained, removing fast time-scale electro- magnetic waves. A. Eulerian Gyrokinetic variational principle We proceed in the reduction of the gyrokinetic variational principle, Eq. (I2TT) . in a very similar way to the Vlasov-Maxwell case (Section HU). The advected parameter is the 6-form f£l = dX A duAdfiAdO B,,f /m = FdXAduAd/uAdO. F is often considered a function for simplicity of notation, but it is understood that operations should be carried out as for a 6-form, e.g., £, V F = V • (FU^) rather than £, V F = U • VF. The connection to the standard distribution function is provided by the Liouville theorem; the advection equation d t (fQ) + £ u (fa) = 0, (23) coupled with Liouville's theorem, Eq. (flTT ). implies d t f + £ v f = 0, (24) which is the standard Vlasov equation. Operating on Eq. (|2~TT ) on the right with the particle-relabeling map, iff' 1 , leads to the reduced 13 Lagrangian l GK (u, ( f> 1 ,A 1 ,F) = L GK ,W~ l ,(f>uA u / dXdudndQF j-A f • t/ x + — fiU e - Hj + -^Jdx (| V0! | 2 - |V x (A + A0| 2 ) , (25) where i/x and Uq are the X and 6 components of the Eulerian fluid velocity U. The unconstrained variations in the Lagrangian frame, 8i]/ lead to constrained variations in the Eulerian frame by defining^S ij(z,t) = 5 l ff(z ,t), (26) or equivalently tj = di/ri/r 1 . Recalling U (z, t) = \j/ (zo, 0, one then calculates drj/dt and 8U, giving Sifr(z ,t)= — 5- — +z ; a ,■ , (27) a? az ; cty ( Zo , = SU (z, + fc j U ^; t \ (28) which is solved for SU = ^ + U • V,?/ - tj • V Z U, (29) giving the variational form stated in Eq. ©. Similarly, using fQ. (z) = / ^o (zo) and 5 (/oQo) = 0, one obtains S(fCl) = -£ n (fCl). (30) Using Eqs. (|29l) and (1301) . we give a basic derivation of the Euler-Poincare equations for a general 14 Lagrangian with an advected volume form, P. dtl = J til = I dt 81 \ 1 81 A 81 j dt^^ J dXdudfidoi^-^-^jT]' 5W 8:U J J *6U' rf + rfPd— 6F , , , d 81 n 81 ,81 \ = / dtl £ v — + FV— , ?7 ' { dt8U U 8U SF (31) giving Eq. © since r\ is arbitrary. The two brackets used are defined as (/u, £) g = £ s J dXdudfidO \i£ between a 1-form density fxdz and a vector field and (f, F^ = J^ S J fP between a function / and a volume form P. Integration by parts is used, with boundary terms dropped, in arriving at the third line. Note that P sometimes includes the volume element (lines 1 and 2) and sometimes does not (lines 3 and 4). A more precise derivation can be found in Refs.Q and 18 It is now simple to write down the equations of motion for U, using SIgk _ e _fj^ SI gk _ mc p 8U X c ' 8Ug e SJgk = 8Jgk = 8U U SU. The derivation is carried out without assumptions about the form of U (e.g., lack of 6 dependence) and the P advection equation [Eq. ([231) 1 is used to cancel time derivatives. We illustrate the general form of the calculation with the 8l/8U' x component of Eq. ©, t — -FAt = — A. — r (FU J ) - -FU J — r dt\c 'I c l dZ J \ I c dZ J c 1 dX e F dX dX \c mc 1 \ + — fiU e - -mu 1 - iuB - H gy \. (33) e 2 The first two terms in Eq. (1331 add to zero due to the advection equation, while the terms involving 15 Ug cancel. Rearranging and expanding the divergence term leads to which gives and -U x xB- mU u b - fiVB - VH gy = 0, (34) U u = -— f ■ (fiVB + VH gv ) (35) mBl v ' B^ c U x = —U x ■ b + — b x (yuVfl + VH gy ) (36) B {1 eB l{ when B 1 - and bx are applied respectively. Similarly, the other 5l GK /6U J equations give U M = 0, (37) 1 dH„ y b-U x = u + -—^, (38) m ou mc mc £ on mc o/u and Eq. (1381) is combined wiht Eq. (1361) to give U x in terms of co-ordinates. The form of Eqs. (1351 - (|39l ) is identical to the standard Lagrangian equations for (x, u,(i, 9} because of the linearity of the Lagrangian in U. Note that various terms have different origins in the Lagrangian and Eulerian derivations; for instance, dX/dt = in the Eulerian derivation (it is just a co-ordinate), while this is not true in the Lagrangian case. Equipped with the solution for U in terms of phase space co-ordinates, the gyrokinetic Vlasov equation is Eq. (1231) . or in co-ordinates, d t F + V- (UP) = 0. (40) Maxwell's equations follow from the standard Euler-Lagrange equations for A and 0, ir = ^r = ' (41) since ^ and A y do not appear in l GK . These lead to the gyrokinetic Maxwell's equations for (pi 16 and Ay, -U 2 «M*) = --^- (42) An 8(p\ (x) J-V x V x A, (x) = 4 VxB - ( 43 ) 4;r oAi (jc) An where *H = Y, s J dXdud/udOFH. IV. THE HAMILTONIAN FORMULATION AND GYROKINETIC POISSON BRACKETS We now perform the generalized Legendre transform of l GK and use the corresponding Hamil- tonian formulation to construct the Poisson brackets for the gyrokinetic system. One complication is the degeneracy in the Lagrangian that arises from the lack of quadratic dependence on U, A and <p. This issue is discussed in detail Refs. LL and gyrokinetic case. 61 and those same arguments apply to the For the moment, we formulate a bracket on the space of plasma densities (see Section IIVBI) and carry out a Legendre transform in U by defining M = SI GK 6U (44) This type of formulation treats the gyrokinetic Poisson- Ampere equations [Eqs. (1421) and (1431 )1 as constraints on the motion of F, rather than dynamical equations in their own right. The gyrokinetic Hamiltonian is defined, as for a standard Legendre transform, as h GK = (M, U\-1 GK (U,P) M-U-F[-A J -U x = ^ J dXdud/udO + ™ nUe-H)}-^ j dx (|V0i | 2 - |V x (A + A0| 2 ) . (45) 17 It is easy to show that with this Hamiltonian [T,&) LP = -(M, 6T SO sm'sm + ( F • V ■ V — ) \ 'SM 6F SM SFlv (46) is a valid Poisson bracket (see Sec. Ill Al and Ref. 1). To evaluate functional derivatives of k [Eq. (|45l) 1. one should obtain the Green's function solutions for <p y and Ay, for instance GK (f ){ (x) = Y J j dX'du'dn'dffK{x\z')F{z), (47) from the gyrokinetic Poisson- Ampere equations, and insert these into h GK , see Refs.l27landl6ll For practical calculation, this is the same as neglecting the electromagnetic part of h GK in the functional derivative. In the same way as the Maxwell- Vlasov system (Sec. Ill A|) . the manifestation of the bracket in Eq. (|46|) is not well defined due to the constraints on M. In the next section, the Dirac theory of constraints (see Appxs. [A]and|B]) is used to reduce Eq. (|46|) to a bracket of the space of densities F. A complete treatment of the geometry of the Poisson- Vlasov system, with the electric field as a constraint, is given in Ref. (6jJ. Many similar ideas will apply to the gyrokinetic system, with complications arising from the nonlocal nature of the theory 21 and larger constraint space (0i and Ai rather than just 0). We reiterate that there are two sets of constraints we consider here; the constraints on M variables, similar to the Maxwell- Vlasov system, and the constraints due to <p\ and Ai, which are the gyrokinetic Poisson- Ampere equations. We first deal with the M constraints, eliminating these variables entirely, then explain how to include <p Y and A l in Section ITVBl A. Gyrokinetic Poisson bracket There are six constraints given by <D / (z) = M / (z)--^- = 0, (48) 8U l (z) 18 with the functional derivatives as listed in Eq. (|32T ). One then forms the constraint matrix Cu (z, z') = {O/ (z) , ®y (z')} with the Poisson bracket of Eq. (g6]) using 5 J ,6{z-z) 6 SS >, 8®i(z) 6Mj(z) — = — AT (z) 5 (z - z) 6 SS , , 6F(z) c <50„(z) sa>, 5F{z) 8F 0, SF(z) -fi6(z- z) 5„< (49) Dropping boundary terms in integrations and inserting the constraint equations (after calculation of the brackets, see Appx. |B]) leads to the very simple form, C 1J {z,z') = F5{z-z')8 ssl x - e -B\ c .' —mb x e -Bl c c - 1 -mb y — W v --b\ c - 1 —mb z SS^ e * mb x mb z e -* e y -<J1£ W e * mc e mc e (50) where all functions are of the z variable and W = R + ^bbVxb. Because of the simple form in z', this matrix is easy to invert according to Eq. ()A2I) giving, u (z, z) = (z - z') <5 M < ' c -b z - c -by ±B\ % e / % - c -b x ^flt £ w z e c - X -B\ m A m -b] mc || — L^t e J e * m -— fl I mc | (51) where again functions are of the z variable, W = b x W and = fi t • W. Of course, this matrix is nothing but the single particle gyrokinetic Poisson matrix 2 as was the case for the Maxwell- Vlasov 19 system. Restricting the functionals Y and to not depend on M (see Appx. |B]) and using d 8T d?SF~ {nFl®j(z)}=F(z)^, (52) the field theory gyrokinetic Poisson bracket is simply, 6F 80\ ,SF 8F) S Here { , } sp is the single particle Poisson bracket structure 2 33 {T,e} DB = {F,\——\ ) . (53) spl v i/.«U = -4 4 -fair- v «r p eBl mBl \ du du + e \ 1 86 8 de) m [dude dude + J. (54) mc \ d/u de de dfi J We note that, although all 8/8M terms are left out above for clarity, the 8/8M u and 8/SM^ terms cancel in a full calculation as expected, so there is no issue with these being undefined (see Appx. |B]). The field theory bracket, Eq. (1531) . is of exactly the form one would expect based on the Poisson- Vlasov bracket 27 and Maxwell- Vlasov bracket— 1 ^ (Eq. (fT3l without 6/6E and 6/6B terms). It is aesthetically pleasing to see this type of structure emerge from the entirely systematic procedure applied above. Since we have not given a proof of the modified Dirac procedure used in the calculation (see Appx. |B]), we directly prove the Jacobi identity for the gyrokinetic bracket [Eq. (l53l) l in Appx. O This shows it is indeed a valid Poisson bracket for this gyrokinetic sys- tem. The reduced Hamiltonian to be used with Eq. (1531) is simply Eq. (|45l) with constraints on M inserted explicitly, h = YjJ dXdud/udeFH - ^- J dx (|V(^| 2 - |V x (A + A0| 2 ) . (55) It is easy to show that d,F = |f, /z} is just the conservative form of the gyrokinetic Vlasov equation, Eq. ii). 20 B. Inclusion of electromagnetic fields The bracket, Eq. (1541 ). does not include electromagnetic field equations, meaning the gyroki- netic Maxwell's equations, Eqs. (1421) and (|43l) , must be specified as separate constraints on the motion to obtain a closed system. Here, we illustrate how to explicitly include the electromagnetic potentials in the bracket for a simplified gyrokinetic system. This procedure also works to extend the simple Poisson-Vlasov bracket 27 61 to include the motion of (p. The general technique is to add a Poisson-Ampere canonical bracket to the gyrokinetic bracket [Eq. (1541 )1 and apply standard Dirac theory to this extended bracket. It is important to recognize that this is only valid because the full constraint matrix would be block diagonal if the reduction were performed in one-step from an original bracket that included electromagnetic and plasma components (i.e., Eq. (l46l) with the addition of canonical brackets in Ay and (py). This condition is satisfied because Ay and <f>\ do not appear in the symplectic structure of the original Lagrangian. For clarity, we use with a simplified electrostatic system in the drift kinetic limit, with H = ecpi + m \6u E \ 2 /2 where 8u E = c(b x V<py) /B. We also assume quasineutrality, which amounts to neglecting the J dx |V0| 2 fSn term in the Lagrangian, and set W to zero 62 . The Hamiltonian for the system is h = / dX(p(X)U(X) Zf * ( m n mc 2 t\ / dXdud/udeF [-u 2 +/uB + e(p + — \V ± <p\ 2 1 , (56) and the unreduced Poisson bracket {r, ®} = (FA —x, —x > + / dX I — — — — — — — J , 57 \ [6F 6F) sp l v J \6cf>6U 6cf) 511) where II = 8118(f) is the variable canonically conjugate to (p. This model is the electrostatic version of the simplified gyrokinetic system in Refs. |57 and 58. Physically, the m\8u E \ 2 /2 = mc 2 /2B 2 |V X 0| 2 term in the Hamiltonian is the polarization drift in the drift kinetic limit 5 ' 58 . V ± indicates a gradient with respect to a co-ordinate system locally perpendicular to the background magnetic field (we are neglecting derivatives of b). Unlike in the previous section, (f> is now con- sidered a separate field in the Hamiltonian, and Poisson's equation should not be used to evaluate functional derivatives. 21 The constraints are 6h vh f ®i = — = ) / dudixdQ eF - mc z V ± ■ ( — V x </> o 2 = n, (58) where II = since 61/ '6<p = 0. $i = is the gyrokinetic Poisson equation; this constraint arises as a secondary Dirac constraint that is necessary to satisfy <t> 2 = 0, see Ref. for more information. Using Eq. (1571 ) the constraint matrix is, Cu(X,X') = c -C (59) where c = ^(Z) 8<P {X') n(X') B 2 (X') i-Y ± 6(X-X') (60) with h = Yjs J dudjidOmF and V ± indicating the derivative is with respect to X' ± . The inverse matrix, chosen to satisfy Eq. (|A2I) . is Cjj(X,X') -c C" 1 (61) where c- l (x,x , ) = --v- ± 1 B 2 (X) h(X) V^6{X-X') (62) The notation V" 1 • [b 2 (X) /h (X) V^ 1 • j is used for clarity to denote the inverse of the operator V ± • [n (X) IB 2 (X) Vjl'L which is similar to the perpendicular Laplacian. While we will not con- sider this here, existence of the inverse of this type of operator (with smooth h/B 2 ) can very likely be proven using standard partial differential equation analysis techniques 63 . The Dirac bracket is 22 constructed using {r,o l( x)} = ^v ± .|^v ± ^ +Nr + —V x . i-^Nr). {r,o 2 (Z)j 6<t>(xy where N r = y / -V • (cfb xV? + £/fltiL ^7 m \ SF m duSF) with the corresponding definition for N & . With Eq. (|A3I) . this leads to dXdudndOF {— , — \ (63) (64) , , B 2 ,8® «5T 501 mc 2 N r + — V ± • (> mc 2 _ JV + — V x • >• (65) Here, { , } sp is the single particle bracket as in Eq. (|54l (with W = 0). Note that in forming Eq. (1651 ), terms involving in the Dirac part of Eq. (IA3I) . canceled with the canonical part of the original bracket, as would be expected. With the reduced Hamiltonian (Eq. (1561 ) without the first term), the bracket can easily be checked to give the Vlasov equation as d t F (X) = yF (X) , h\ ■ Noticing that the d/du term in d t F [see Eq. (1401) 1 integrates to zero, we see that f dudfxdO ed t F = -Nh- This is used in d t <p(X) = {<f>(X),h} DB to show y J dudfidO i^ed t F - mc 2 V ± B 2 dt V 0, (66) which is just the time derivative of Poisson's equation for this gyrokinetic model. Using the procedure presented above there should be no particular obstacle to the construction of brackets for more complex gyrokinetic theories. For instance, one could include finite Larmor radius effects or magnetic fluctuations 64 . However, considering the complexity of the bracket for even a very simple gyrokinetic model, such brackets are unlikely to be of much practical use. 23 V. NUMERICAL APPLICATIONS In a general sense, one of the main motivations for this work is to possibly help address the increasing need to develop new algorithms with the long-time conservation properties necessary to help improve the physics fidelity of simulation results as we move towards exascale computing and beyond. Of specific interest in this paper is the desire to explore the possibility of applying the Eulerian variational methods developed here in a discrete context to design continuum geometric integrators for gyrokinetic systems. To elaborate on this idea, here we give a simple explanation of some geometric discretization methods based on recent work in numerical fluid dynamics. The methods described here are just examples from a large array of literature on the subject. Some 65 69 In addition we remark on how Euler- other techniques can be found in, for instance Refs. Poincare models can be used to formulate sub-grid models for turbulence simulation and some of the challenges associated with extending these ideas to gyrokinetic turbulence. a. Lagrangian side: discrete Euler-Poincare equations Conceptually, an obvious way to design a geometric integrator for an Euler-Poincare system is to directly discretize the Euler- Poincare variational principle. If one can design discrete variations of the correct form, the entire integrator can be constructed directly from the variational principle as for a standard variational integrator. This approach has recently been successfully applied to develop an integrator for the Euler fluid equations 44 and more complex fluids, including magnetohydrodynamics (MHD) 45 . The utility of such an approach is illustrated by the very nice properties of these schemes. For instance, the MHD scheme 45 exactly preserves V B = and the cross helicity J dx V'B. As one consequence of this, there is almost no artificial magnetic reconnection. The symplectic nature of the scheme also leads to other very good long time conservation properties. The first requirement in constructing an Euler-Poincare integrator is a finite dimensional ap- proximation to the diffeomorphism Lie group. In the case of fluids or MHD, the group is that of volume preserving diffeomorphisms and a matrix Lie group is constructed to satisfy analogous properties to the infinite dimensional group. For Vlasov-Poisson, Vlasov-Maxwell or a gyrokinetic system, the group is that of symplectomorphisms. Thus, for a discretization, a different matrix Lie group than the fluid case should be used, with properties designed to mimic those of the infinite dimensional symplectomorphism group. Using this group one can find the Lie algebra, which will give the form of the space of discrete vector fields (just the Eulerian phase space fluid veloc- ities). Group operations can then be constructed as matrix multiplications as for a standard finite 24 dimensional Lie group, and advected parameters included through the use of discrete exterior cal- culus. One would then use the discrete Euler-Poincare theorem 10 , which gives discrete update equations [in analogy with Eq. (|3TT) 1 from a discrete reduced Lagrangian. An algorithm of this form can be shown to be symplectic and have similar conservation properties (arising from vari- ants of Noether's theorem) to the continuous system. The final update equations obtained from this method are not as complex as one might expect and would not preclude incorporation into large scale codes. Obviously there are several unanswered questions regarding the application of this method to kinetic plasma systems. First, one must discretize the symplectomorphism group, which may not be trivial. In two phase space dimensions the group is the same as the group of vol- ume preserving diffeomorphisms; however, in higher dimensions the symplectomorphisms form a more restricted class of transformations. The lack of a finite boundary in velocity space may also present issues relating to the discretization of the symplectomorphisms. The degeneracy of the system is another aspect which differs from the fluid system, and the consequences of this in the discrete setting would have to be carefully considered. Finally, for a gyrokinetic system, it would be necessary to remove the 6 dimension in some way. This could potentially be done either in the continuous setting or after discrete equations have been obtained. b. Hamiltonian side: Poisson bracket discretization Another way to form a discrete Hamil- tonian system is to directly discretize the Poisson bracket. The general idea is simple, one finds a discrete Hamiltonian functional and discrete bracket that are finite dimensional approximations to the continuous versions. In this way, one discretizes (in phase space) via the method of lines, and reduces the infinite dimensional system to an approximate finite dimensional one. Any sym- plectic temporal discretization can then be used to ensure the system is discretely Hamiltonian 69 . The difficultly arises in ensuring a correct Hamiltonian discretization of the bracket. This requires antisymmetry and the Jacobi identity to be satisfied, and such a bracket can be very difficult to find in practice. For instance, for the Euler fluid equations, the non-canonical structure complicates matters and a discrete bracket has been found only for simplified cases^-. An obvious place to start in this endeavor would be the Vlasov-Poisson system, as the structure is much more simple. Generalizations to gyrokinetic systems could then potentially be achieved through Nambu bracket formulations 57 . c. Alpha models and large-eddy simulation Much work has been done in the last decade in the fluids community on so-called alpha models. The general idea is to regularize the fluid equa- tions (Navier-Stokes or MHD) at the level of the Euler-Poincare variational principle, by adding 25 terms into the Lagrangian that include gradients of the fluid velocity. These terms penalize the formation of small scale structures, and can thus be used as a large eddy simulation (LES) model, causing turbulence to dissipate at larger scales^. These methods have been shown to have some significant advantages over more traditional LES methods (for instance those based on hyperdif- fusion) especially for simulation of MHD turbulence3Z2. As gyrokinetic turbulence simulation becomes a more mature subject, it is interesting to enquire whether similar alpha models could be formulated for gyrokinetic large eddy simulation. In fact, alpha models can be derived from a standard fluid variational principle, by averaging over small scale fluctuations that are assumed to be advected by the larger scale flow^^. Ap- proaching the gyrokinetic variational principle in a similar way leads to the addition of extra, regu Ref. tromagnetic gauge invariance, we were led to the regularized Lagrangian I = l GK + l a , where arized terms into the gyrokinetic Lagrangian. For instance, following the general ideas in 5l[ averaging over perpendicular X-space fluctuations of scale length a, and ensuring elec- l a = a 2 J dXdudfideP • V x U x - V 2 ± h) . (67) While this Lagrangian gives well-defined equations of motion, there is a fundamental problem in that it destroys some of the degeneracy in the original system. As a consequence of this, the equations of motion involve solving spatial PDEs for U, which would significantly increase com- putation times, defeating the purpose of an LES. It is not yet clear if it is possible to design a regularization of this type for the gyrokinetic system that retains the redundancy of the U fields and allows one to write down a standard Vlasov equation. We note that gyrokinetic LES has been explored and implemented recently on the GENE code, by adding hyperdiffusive terms in the perpendicular co-ordinate s 49174 . VI. CONCLUDING REMARKS In this article we have applied the Euler-Poincare formalism to derive a new gyrokinetic action principle in Eulerian co-ordinates. We start with a single-particle Poincare-Cartan 1-form, using the theory of Ref. |21| to systematically construct a gyrokinetic field theory action in Lagrangian co-ordinates. The fundamental idea is then to reduce this action using symmetry under the the particle-relabeling map, (jc, v) = i/r(jc ,v ), which takes particles with initial position (x ,v ) to 26 their current location (x, v) 1 ' 18 ' 19 . This process leads to an action functional formulated in terms of the Eulerian phase space fluid velocity, U, and the advected plasma density, F, as well as the standard electromagnetic potentials. In the course of reduction, the arbitrary variations of the La- grangian fields (used to derive equations of motion) lead to constrained variations of the Eulerian fields, U and F. Because of this, field motion is governed by the Euler-Poincare equations rather than the standard Euler- Lagrange equations. Explicit calculation of the Euler-Poincare equations for a standard gyrokinetic single particle Lagrangian is shown to give the gyrokinetic Vlasov equation. Since the space of electromagnetic potentials is not altered by iff (xq, Vq), the gyrokinetic Poisson-Ampere equations arise from the standard Euler-Lagrange equations for the perturbed potentials. Using the methodology set out in Ref. 1 we then perform a Legendre transform to derive the Hamiltonian form of the gyrokinetic system. The principal difficulty is the strong degeneracy, which is related to the linearity and lack of time derivatives for certain function variables in the action principle. Physically, this arises from the fact that the plasma distribution function encodes the information about particle phase space trajectories. The degeneracy leads to a Poisson bracket in terms of too many variables; namely, a series of constrained momentum variables canonically conjugate to U as well as the distribution function F. To reduce the bracket into a well defined form we use a modified version of the Dirac theory of constraints, which is a systematic way to project a Poisson bracket onto a constraint submanifold when momentum variables are constrained. The modified Dirac procedure (see Appx. EJ can be a significant simplification over standard Dirac theory for certain types of systems. This process leads to an infinite dimensional gyrokinetic Pois- son bracket, which takes a natural form based on the single particle bracket. We also demonstrate how this procedure leads to the full, electromagnetic Vlasov-Maxwell bracket 2 2^^. Since the electromagnetic equations in the gyrokinetic system are really constraints on the motion, we chose to include these in the bracket via a second application of the Dirac theory of constraints. The general method is expounded through construction of the bracket for a simplified electrostatic model with no finite Larmor radius effects. Although the brackets obtained by such an approach are probably to complicated to be of much practical use, it makes for an interesting application of Dirac theory. As a final comment, we note that there is an emerging (and very likely increasing) demand for a new class of algorithms capable of dealing with the demands of powerful modern super- computers - moving on the path to the exascale and beyond 75 . Approaches based on geometric 27 integration may prove highly relevant as increasingly ambitious simulations of larger and more complex systems are undertaken. ACKNOWLEDGEMENTS We extend our thanks to Dr. John Krommes and Joshua Burby for enlightening discussion. This research is supported by U.S. DOE (DE-AC02-09CH1 1466). CC acknowledges financial support from the Agence Nationale de la Recherche and from the European Community under the contract of Association between EURATOM, CEA, and the French Research Federation for fusion studies. Appendix A: Dirac Constraints The Dirac theory of constraints or Dirac bracket is used to build Poisson brackets for Hamil- tonian systems with constraints. The original purpose of the theory was to construct quantizable Poisson brackets starting with a degenerate Lagrangian, i.e., a Lagrangian where the momenta are not independent functions of velocities. The theory applies equally well to a bracket that is already in non-canonical form, a realization that can be very useful in the construction of field theoretic brackets 36 - 53 . For example, in Ref.[35|, the non-canonical magnetohydrodynamic bracket is reduced to incorporate the incompressibility constraint. We give a very brief overview of the theory here for the convenience of the reader. More complete treatments can be found in Refs. 76. 34H36J, 53 and We consider an infinite dimensional Hamiltonian system, described by field variables Xi(z) (for i = 1, . . .,«), with Poisson bracket { , }, Hamiltonian H, and m local constraint functions ®\(z), 2 (z), . . . , 0/v(z) = 0. The constraint matrix, C i} (z,z') = {o i (z),0j(z')}, (Al) and its inverse, defined using dz'dj (z, z) C]l (z', z") = 5 ik 8 (z - z") , (A2) 28 are used to form the Dirac bracket, {r,0} DB = {r,©} - J dzdz' {r,^(z)}C u l (z,z'){®j(z'),®} (A3) By construction this bracket satisfies the Jacobi identity^^^. Geometrically, the constraints force motion to lie on a constraint submanifold, which inherits the Dirac bracket from the Poisson bracket on the original manifold 76 . More precisely, the constraints O are Casimir invariants of the Dirac bracket (|A3I) . i.e., {r, 0,} D b = f° r an y functional Y and for all j = 1, . . . , m. In the case where the matrix C is not invertible, Dirac theory suggests the use of one or more secondary constraints, which must be included and the constraint matrix C recalculated. See Refs. 35 and 76 for more information. Appendix B: Modified Dirac procedure In this Appendix we provide some justification of the modified Dirac procedure used in the calculation of the Vlasov-Maxwell and gyrokinetic Poisson brackets [Eqs. (|T3l) and (l53ll. For clarity we restrict ourselves to the finite dimensional case, even though the technique is applied to infinite dimensional systems in the present manuscript. The purpose of this modified procedure is to simplify the computation of the Dirac bracket and reduce the dimensionality of the constrained system. The modified Dirac procedure involves the following steps: 1. Compute the constraint matrix C and simplify its expression by applying the constraints. The resulting matrix C is weakly equal to C. 2. Construct the Dirac bracket using Eq. (IA3I ) with C instead of C. 3. Choose a set of m variables, denoted z, to be eliminated from the bracket (IA31 ) using the constraint equations, so that it can be completely rewritten in terms of the remaining n - m variables, denoted z. Substitute these into the bracket and drop all partial derivatives with respect to the z. This new bracket (of reduced dimensionality) is the reduced Dirac bracket. Since the above procedure departs from the standard one for the computation of the Dirac bracket, 29 the reduced bracket is not a priori a Poisson bracket. Below we justify this procedure by providing the condition under which the reduced Dirac bracket is a Poisson bracket. Regarding Step 1, assuming C is analytic in the constraint variables, the invertibility of the matrix C is sufficient to ensure that the reduced Dirac bracket exists, in particular, that it satisfies the Jacobi identity. This follows from the fact that if det C = det [C (O = 0)] is non-zero, we can find a small open neighborhood of O = such that det C ^ by continuity, implying C is invertible in this region. Regarding Step 3, we write the Dirac bracket in the form oz oz where is the Poisson matrix 24 satisfying the Jacobi identity, = ^ Jpit" + ft J™ + . (B 1 ) As explained above, we single out a subset of the z variables to be eliminated by the reduction procedure (denoted z) and label the remaining variables z. We rewrite the m constraints in the form 0;(z) = Zj ~ iPjil). This should be possible at least locally under the hypothesis of the implicit function theorem. Note that all the constraints considered in this manuscript are already of this form. Next we consider the general coordinate change z = (z, z) i-> w), where tD are the constraints. This change of variables is invertible, and its inverse is given by z,j = O y - + <pj(w) and z = w. Since the constraints are Casimir invariants of the Dirac bracket, the transformed Jacobi matrix must be in the form '0 "D = (B2) D J where the number of rows and columns of zeros is the same as the number of constraints. Since I D satisfies the Jacobi identity, it is straightforward to see from Eq. (|B1I) that Y D must also. Setting O to and applying the coordinate change whz leads to a Poisson bracket in z that satisfies the Jacobi identity. The set of functions depending on z constitutes a Poisson algebra with the Poisson 30 bracket given by [f > 8} ' = di- 3 °di- This bracket is obtained by simply dropping the d/dz terms in the Dirac bracket. For con- sistency we require this to be equal to the operation of the reduced Dirac bracket {f r ,g r }' DB , a condition which is obviously ensured if we drop d/dz and use z = z w) to eliminate z. Note that this is equivalent to removing the rows and columns of J D corresponding to the z variables. Of course, one can reach the same point by the variable change z i-> w) as explained above, illustrating that dropping z partial derivatives will lead to a bracket satisfying the Jacobi identity. Note that if a constraint is simply a co-ordinate O, = zj, all d/dzj will be automatically elimi- nated from the Dirac bracket by the standard Dirac procedure. This behavior is seen in both the finite dimensional example below as well as the infinite dimensional brackets in the main body of the paper. For example, the functional derivative S/6M V cancels in the calculation leading to the Vlasov-Maxwell bracket, Eq. (PT3l . We now give an example of the procedure for a well known finite dimensional system, the Poisson bracket for the Lorentz force. Using a 12-dimensional extended phase space (x, v,p x ,p v ) and the modified Dirac theory above, both the (jc, v) bracket and the canonical bracket are easily derived in one step. Start with the Lagrangian for a particle in an electromagnetic field (using unit charge and mass for clarity), 1 7 L = (A + v) ■ x - -v - <f> (B3) and apply the Dirac theory of constraints to the canonical Poisson bracket in extended phase space df dg df dg dxi dp xi dvi dp vi where v. v. means switch f to g and vice versa. There are 6 constraints on the canonical momenta, ®xi - Pxi ~ (Aj + Vi) = 0, <D V! = p vi = obtained from dL/dii = p z \. This gives the constraint matrix C = d(A + v), the symplectic form in (x, v) space. The Dirac bracket follows from the computation of C 1 : {fg} = £^J^_ VV +f^_ vv DB dxj dp xi ' dxj dv{ df dgdAi „ df dg + ir-irir 1 - v - v - + B -ir x ir- ^ dp X j dVj dxj dv dv 31 Note that, as explained above, the terms involving d/dp V j cancel automatically. We can now form the reduced Dirac bracket in the variables of our choice (except p v ) by simply removing partial derivatives and substituting in constraints. Here there is no need to apply the constraints on C before its inversion, so Step 1 has been skipped. For instance, removing p x trivially leads to the standard (x, v) bracket [f(x,v),g(x,v)} DB = —— - v.v. + B ■ — x—, (B5) OXj OVi ov ov while removing v leads to the (x, p x ) canonical bracket, {f(x,p x ),g(x,p x )} DB = ^ d -v.v. (B6) One could even form a (v,p x ) bracket if desired, although the dAj/dxj term would have to be changed into (v,p x ) co-ordinates using the constraint equations. While these finite dimensional brackets can be straightforwardly derived from the Lagrangian [Eq. (|B3I )1 through other means, derivation of infinite dimensional brackets can be more challeng- ing. So long as suitable care is taken, the method outlined above can be very useful in deriving Poisson brackets from certain types of field theoretic Lagrangians, as is carried out in the main body of this work for the Vlasov-Maxwell and gyrokinetic brackets. Appendix C: Jacobi identity for the gyrokinetic Poisson bracket We provide a direct proof of the Jacobi identity, {r,{©,A}}+ 0=0, (ci) for the gyrokinetic bracket presented above, (6T_ 6®\ {r,0}= / dzF{—, — \ . (C2) \^F , 6F) sp Here O denotes the permutation of Y, 0, and A through the other two possibilities. The bracket, Eq. (IC2I) . is in the form / 6F \SF 32 £ S 4), (C3) 5F where is an anti-self adjoint operator with dependence on P. For operators of this type (see Ref. 27) the second functional derivatives cancel in the calculation of Eq. (IC1I) . This implies that the only part of 6{&, A} /6F that contributes in the Jacobi identity is IST/dF, dQ/dF\ and l I sp Eq. (IC1I) becomes r , (sr (se SA] } , n H{w{wvll +Um0 - <C4) Thus, the Jacobi identity of the functional Poisson bracket follows directly from that for the single particle bracket. REFERENCES ! H. Cendra, D. D. Holm, M. J. W. Hoyle, and J. E. Marsden, Journal of Mathematical Physics 39,3138 (1998). 2 A. Brizard and T. Hahm, Reviews of Modern Physics 79, 421 (2007). 3 J. Krommes, Annual Review of Fluid Mechanics 44, 175 (2012). 4 Z. Lin, T. Hahm, W. Lee, W. Tang, and R. White, Science 281, 1835 (1998). 5 D. Dubin, J. Krommes, C. Oberman, and W. Lee, Physics of Fluids 26, 3524 (1983). 6 R. G. Littlejohn, Journal of Mathematical Physics 23, 742 (1982). 7 T. S. Hahm, W. W. Lee, and A. Brizard, [Physics of Fluids 31, 1940 (1988) 8 H. Sugama, Physics of Plasmas 7, 466 (2000). 9 A. J. Brizard, Physics of Plasmas 7, 4816 (2000). 10 P. Similon, Physics Letters A 112, 33 (1985). n B. M. Boghosian, Covariant Lagrangian Methods of Relativistic Plasma Theory, Ph.D. thesis, University of California, Davis (1987). 12 D. Pfirsch and P. J. Morrison, Physical Review A 32, 1714 (1985). 13 D. Pfirsch and P. J. Morrison, Physics of Fluids B 3, 271 (1991). 14 H. Ye and P. Morrison, Physics of Fluids B 4, 771 (1992). 15 A. Brizard, Physical Review Letters 84, 5768 (2000). 16 F. Low, Proceedings of the Royal Society of London A 248, 282 (1958). 17 T. Fla, Physics of Plasmas 1, 2409 (1994). 18 D. Holm, T. Schmah, and C. Stoica, Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions, Oxford Texts in Applied and Engineering Mathematics (Oxford University Press, 33 USA, 2009). 19 D. D. Holm, J. E. Marsden, and T. S. Ratiu, Advances in Mathematics 137, 1 (1998). 2() W. A. Newcomb, Nuclear Fusion: Supplement, part 2 , 451 (1962). 21 H. Qin, R. H. Cohen, W. M. Nevins, and X. Q. Xu, Physics of Plasmas 14, 0561 10 (2007). 22 J. E. Marsden and A. Weinstein, Physica D 4, 394 (1982). 23 P. J. Morrison, [Physics of Plasmas 20, 012104 (2013)| 24 C. Chandre, L. D. Guillebon, A. Back, E. Tassi, and P. Morrison, (2012), arXiv: 1205.2347 [math-ph]] 25 P. Morrison, Physical Letters A 80A, 383 (1980). 26 P. Morrison and J. Greene, Physical Review Letters 45, 790 (1980). 27 P. J. Morrison, in A merican Institute of Physics Conference Series} American Institute of Physics Conference Series, Vol. 88, edited by M. Tabor and Y. M. Treve (1982) pp. 13-46. 28 P. Morrison, M. Vittot, and L. D. Guillebon, (2012), |arXTv: 1212.3 007 [ph ysics.plasm -ph] . 29 D. Pfirsch and D. Correa-Restrepo, Journal of Plasma Physics 70, 719 (2004). 3() D. Correa-Restrepo and D. Pfirsch, Journal of Plasma Physics 71, 503 (2005). 31 D. Correa-Restrepo and D. Pfirsch, Journal of plasma physics 70, 199 (2004). 32 D. Correa-Restrepo and D. Pfirsch, Journal of Plasma Physics 70, 757 (2004). 33 D. Correa-Restrepo and D. Pfirsch, Journal of Plasma Physics 71, 1 (2005). 34 P. M. Dirac, Canadian Journal of Mathematics 2, 129 (1950). 35 C. Chandre, P. Morrison, and E. Tassi, Physics Letters A 376, 737 (201 1). 36 C. Chandre, E. Tassi, and P. J. Morrison, Physics of Plasmas 17, 042307 (2010). 37 J. E. Marsden and M. West, Acta Numerica 10, 357 (2001). 38 H. Qin and X. Guan, |Phys. Rev. Lett. 100, 035006 (2008)1 39 J. Li, H. Qin, Z. Pu, L. Xie, and S. Fu, [Physics of Plasmas 18, 052902 (20 lT) 40 J. Squire, H. Qin, and W. M. Tang, |Physics of Plasmas 19, 052501 (20T2)| 41 A. Lew, J. E. Marsden, M. Ortiz, and M. West, Arch. Rational Mechanics and Analysis 167, 85 (2003). 42 J. Marsden, G. Patrick, and S. Shkoller, Communications in Mathematical Physics 199, 351 (1998). 43 A. Stern, Y. Tong, M. Desbrun, and J. E. Marsden, (2007), |arXTv:0707.4470v3 [math.NAfl 44 D. Pavlov, P. Mullen, Y. Tong, and E. Kanso, Physica D 240, 443 (2009). 45 E. Gawlik, P. Mullen, D. Pavlov, J. Marsden, and M. Desbrun, Physica D 240, 1724 (2011). 34 46 J. Squire, H. Qin, and W. M. Tang, Physics of Plasmas 19, 084501 (2012). 47 D. D. Holm, J. E. Marsden, T. Ratiu, and A. Weinstein, Physics Reports 123, 116 (1985). 48 M. Kruskal and C. Oberman, Physics of Fluids 1, 275 (1958). 49 P. Morel, A. Navarro, M. Albrecht-Marc, D. Carati, F. Merz, T. Gorier, and F. Jenko, Physics of Plasmas 18, 072301 (2011). 50 J. P. Graham, D. Holm, P. Mininni, and A. Pouquet, Journal of Scientific Computing 49, 21 (2011). 51 H. Bhat, R. Fetecau, J. Marsden, K. Mohseni, and M. West, Multiscale Modelling and Simula- tion 3, 818 (2005). 52 D. D. Holm, Chaos 12, 518 (2002). 53 P. Morrison, N. R. Lebovitz, and J. A. Biello, [Annals of Physics 324, 1747 (2009)| 54 D. Holm and C. Tronci, Communications in Mathematical Sciences 10, 191 (2012). 55 A. J. Brizard, |Journal of Plasma Physics 71, 225 (2005)] 56 C. Tronci, [Plasma Physics and Controlled Fusion 55, 035001 (2013")j 57 B. Scott, A. Kendl, and T. Ribeiro, Contribributions to Plasma Physics 50, 228 (2010). 58 B. Scott and J. Smirnov, Physics of Plasmas 17, 1 12302 (2010). 59 If space-time is curved, the d/dt component of r could have dependence on the phase-space co-ordinates. When this is true, one can not work in 6-D phase space with a separate treatment of time. 6() C. Chandre, A. Brizard, and E. Tassi, (2012), |arXiv: 121 1.0850vl [nlin.CD] 61 H. Giimral, J. Math. Phys. 51, 083501 (2010). 62 With W = the gyrokinetic Lagrangian is no longer gyrogauge invariant, an issue of little practical importance for formulating such a model. 63 L. Bers, F. John, and M. Schechter, in Proceedings of the summer seminar in Boulder Colorado, Lectures in applied mathematics, Vol. 3, edited by A. Householder (1957). 64 Including full magnetic fluctuations would require careful treatment of the electromagnetic gauge. This could be handled using first class constraints for unphysical degrees of freedom^ 4 - or by fixing the gauge with a Lagrange multiplier, A, and treating A as a separate variable 8 . 65 C. J. Cotter, D. D. Holm, and P. E. Hydon, Proceedings of the Royal Society of London: A 463, 2671 (2007). 66 S. Elcott, Y. Tong, E. Kanso, P. Schroder, and M. Desbrun, |ACM Transactions on Graphics 26 (2007), 10. 1145/11 89762. 11897661 35 6 7 P. Mullen, K. Crane, D. Pavlov, Y. Tong, and M. Desbrun, ACM Transactions on Graphics 28, 38:1 (2009). 68 X. Zhang, D. Schmidt, and B. Perot, |Journal of Computational Physics 175, 764 (2002)| 69 T. Bridges and S. Reich, Journal of Physics A 39, 5287 (2006). 70 N. Bou-Rabee and J. Marsden, Foundations of Computational Mathematics 9, 197 (2009). 71 R. McLachlan, Physical Review Letters 71, 3043 (1993). 72 S. Chen, C. Foias, D. Holm, E. Olson, E. Titi, and S. Wynne, Physical Review Letters 81, 5338 (1998). 73 S. Chen, D. Holm, L. Margolin, and R. Zhang, Physica D 133, 66 (1999). 74 P. Morel, A. Banon Navarro, M. Albrecht-Marc, D. Carati, F. Merz, T. Gorier, and F. Jenko, Physics of Plasmas 19, 012311 (2012) 75 R. Rosner et al. (DoE Advanced Scientific Computing Advisory Committee Report), "The opportunitie s and ch allenge s of exascale com puting," (2010). 76 J. Marsden and T. Ratiu, Introduction to Mechanics and Symmetry: A Basic Exposition of Clas- sical Mechanical Systems, Texts in Applied Mathematics (Springer, 1999). 36